Question 2: Distribution of Recovery Time from A
Surgery
Time to recovery (in days) after a specific knee surgery procedure.
This follows a typical log-logistic pattern in medical
survival/recovery analysis:
8.23, 12.74, 14.83, 16.61, 18.16, 19.55, 20.80, 21.94, 23.00, 23.98, 24.89, 25.75, 26.56,
27.34, 28.08, 28.79, 29.48, 30.15, 30.81, 31.45, 32.08, 32.70, 33.31, 33.92, 34.53, 35.13,
35.73, 36.33, 36.93, 37.53, 38.14, 38.75, 39.37, 40.00, 40.64, 41.29, 41.95, 42.63, 43.33,
44.05, 44.79, 45.56, 46.36, 47.20, 48.08, 49.02, 50.03, 51.12, 52.32, 53.65
Based on the above data to perform the following analysis.
- Using method of moment estimation to estimate \(\alpha\) and \(\beta\), denoted by \(\hat{\alpha}\) and \(\hat{\beta}\), respectively.
Moment generating function for a log logistic distribution. \[
\mu_k = E[X^k] = \alpha^k B\left(1+\frac{k}{\beta},\,
1-\frac{k}{\beta}\right)
\] The first two moments are \[
\mu_1 = E[X] = \alpha B\left(1+\frac{1}{\beta},\,
1-\frac{1}{\beta}\right)
\] \[
\mu_2 = E[X^2] = \alpha^2 B\left(1+\frac{2}{\beta},\,
1-\frac{2}{\beta}\right)
\] Beta Function Identity \[
B(a,b)=\frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}
\] Simplifying both moments \[
\mu_1
=
\alpha
\Gamma\left(1+\frac{1}{\beta}\right)
\Gamma\left(1-\frac{1}{\beta}\right)
\] \[
\mu_2
=
\alpha^2
\Gamma\left(1+\frac{2}{\beta}\right)
\Gamma\left(1-\frac{2}{\beta}\right)
\] The first two sample moments are \[
m_1 = \frac{1}{n}\sum_{i=1}^n X_i = 34.1922
\] \[
m_2 = \frac{1}{n}\sum_{i=1}^n X_i^2 = 1288.845
\] Substitute theoretical sample moments for moments since
they’re equal \[
34.1922
=
\alpha
\Gamma\left(1+\frac{1}{\beta}\right)
\Gamma\left(1-\frac{1}{\beta}\right)
\]
\[
1288.845
=
\alpha^2
\Gamma\left(1+\frac{2}{\beta}\right)
\Gamma\left(1-\frac{2}{\beta}\right)
\] Solve first equation for \(\alpha\) \[
\alpha
=
\frac{34.1922}
{\Gamma\left(1+\frac{1}{\beta}\right)
\Gamma\left(1-\frac{1}{\beta}\right)}
\]
Plug \(\alpha\) into second equation
\[
1288.845
=
\left(
\frac{34.1922}
{\Gamma\left(1+\frac{1}{\beta}\right)
\Gamma\left(1-\frac{1}{\beta}\right)}
\right)^2
\Gamma\left(1+\frac{2}{\beta}\right)
\Gamma\left(1-\frac{2}{\beta}\right)
\]
Divide both sides by \({m_1^2}\)
(replacing theoretical sample moments with \({m_1}\) and \({m_2}\)) \[
\frac{m_2}{m_1^2}
=
\frac{
\Gamma\left(1+\frac{2}{\beta}\right)
\Gamma\left(1-\frac{2}{\beta}\right)
}{
\left[
\Gamma\left(1+\frac{1}{\beta}\right)
\Gamma\left(1-\frac{1}{\beta}\right)
\right]^2
}
\]
Define \[
g(\beta)
=
\frac{
\Gamma\left(1+\frac{2}{\beta}\right)
\Gamma\left(1-\frac{2}{\beta}\right)
}{
\left[
\Gamma\left(1+\frac{1}{\beta}\right)
\Gamma\left(1-\frac{1}{\beta}\right)
\right]^2
}
-
\frac{m_2}{m_1^2}
\] Solve \(g(\beta)=0\) to
obtain \(\hat{\beta}\), and then
substitute into \[
\hat{\alpha}
=
\frac{m_1}
{\Gamma\left(1+\frac{1}{\hat{\beta}}\right)
\Gamma\left(1-\frac{1}{\hat{\beta}}\right)}
\]
Code to solve \(g(\beta)=0\) and
obtain both \(\hat{\beta}\) and \(\hat{\alpha}\).
x <- c(8.23, 12.74, 14.83, 16.61, 18.16, 19.55, 20.80, 21.94, 23.00, 23.98, 24.89, 25.75, 26.56,
27.34, 28.08, 28.79, 29.48, 30.15, 30.81, 31.45, 32.08, 32.70, 33.31, 33.92, 34.53, 35.13,
35.73, 36.33, 36.93, 37.53, 38.14, 38.75, 39.37, 40.00, 40.64, 41.29, 41.95, 42.63, 43.33,
44.05, 44.79, 45.56, 46.36, 47.20, 48.08, 49.02, 50.03, 51.12, 52.32, 53.65)
#1st and 2nd moments
m1 = mean(x)
m2 = mean(x^2)
#root finding equation
gb <- function(B){
( (gamma(1 + 2/B) * gamma(1 - 2/B)) / ((gamma(1 + 1/B) * gamma(1 - 1/B))^2) ) - (m2 / m1^2)
}
# since beta > 2
x_vals <- seq(2.01, 10, length.out = 200)
dframe <- data.frame(x = x_vals, y = gb(x_vals))
gk.plot <- ggplot(dframe, aes(x = x, y = y)) +
geom_line(color = "steelblue", size = 1) +
geom_hline(yintercept = 0, linetype = "dashed", alpha = 0.5) + # x-axis
geom_vline(xintercept = 0, linetype = "dashed", alpha = 0.5) + # y-axis
labs(title = "The curve of function g(Beta)",
x = "Beta",
y = "g(Beta)") +
theme(plot.title = element_text(hjust = 0.5),
plot.margin = margin(t = 35, r = 20, b = 30, l = 30, unit = "pt"))
ggplotly(gk.plot)
#find the root (beta value when equation = 0)
B <- uniroot(gb, interval = c(5, 7))$root
#use the beta to find the alpha
alpha <- m1/ (gamma(1+1/B) * gamma(1-1/B))
pander(cbind(beta = B, alpha = alpha))
- Since the moment estimates \(\hat{\alpha}\) and \(\hat{\beta}\) are random, construct
bootstrap sampling distributions for each. To visualize these
distributions, plot separate bootstrap histograms for \(\hat{\alpha}\) and \(\hat{\beta}\). hen, overlay a smooth
density curve on each histogram using Gaussian kernel density
estimation. Finally, describe the patterns of these density curves.
bootstrap <- function (data, statistics, B){
n <- length(data)
stats <- matrix(nrow = B, ncol = 2)
for (i in 1:B){
new_data <- sample(data, n, replace = TRUE)
stats[i,] <- statistics(new_data)
}
return(stats)
}
statistics <- function (data){
m1 = mean(data)
m2 = mean(data^2)
gb <- function(B){
( (gamma(1 + 2/B) * gamma(1 - 2/B)) / ((gamma(1 + 1/B) * gamma(1 - 1/B))^2) ) - (m2 / m1^2)}
x_vals <- seq(2.01, 200, length.out = 100)
dframe <- data.frame(x = x_vals, y = gb(x_vals))
B <- uniroot(gb, interval = c(3, 15))$root
alpha <- m1/ (gamma(1+1/B) * gamma(1-1/B))
c(B, alpha)
}
stats <- bootstrap(x, statistics, 5000)
beta_boot <- stats[, 1]
alpha_boot <- stats[, 2]
hist(beta_boot,
breaks = 20,
probability = TRUE,
xlab = expression(hat(beta)),
main = "Bootstrap Sampling Distribution\n of Beta Hat",
cex.main = 0.9,
col = "lightgray",
border = "white")
lines(density(beta_boot), col = "navy", lwd = 2)

The bootstrap sampling distribution of \(\hat{\beta}\) is roughly bell-shaped with a
little bit of skewness to the right. It’s centered close to the original
MoM estimate, which suggests the estimator is not heavily biased. The
spread of the histogram shows how much \(\hat{\beta}\) would vary if we collected
new samples from the population. Since the distribution isn’t extremely
wide or irregular, it suggests that the MoM estimator for \({\beta}\) is fairly stable for this
sample.
hist(alpha_boot,
breaks = 20,
probability = TRUE,
xlab = expression(hat(alpha)),
main = "Bootstrap Sampling Distribution\n of Alpha Hat",
cex.main = 0.9,
col = "lightgray",
border = "white")
lines(density(alpha_boot), col = "darkred", lwd = 2)

The bootstrap sampling distribution of \(\hat{\alpha}\) is approximately symmetric
and bell-shaped. It’s centered pretty close to the original MoM
estimate, which suggests the estimator is not heavily biased. The spread
of the histogram shows how much \(\hat{\alpha}\) would vary if we collected
new samples from the population. Since the distribution isn’t extremely
wide or irregular, it suggests that the MoM estimator for \({\alpha}\) is fairly stable for this
sample.
LS0tCnRpdGxlOiAiQXNzaWdubWVudCAzOiBNZXRob2RzIG9mIE1vbWVudCBFc3RpbWF0aW9uIgphdXRob3I6ICJDaGFybGllIE1vcmdhbiIKZGF0ZTogIkR1ZTogMi8yMy8yNiIKb3V0cHV0OgogIGh0bWxfZG9jdW1lbnQ6IAogICAgdG9jOiB5ZXMKICAgIHRvY19kZXB0aDogNAogICAgdG9jX2Zsb2F0OiB5ZXMKICAgIG51bWJlcl9zZWN0aW9uczogbm8KICAgIHRvY19jb2xsYXBzZWQ6IHllcwogICAgY29kZV9mb2xkaW5nOiBoaWRlCiAgICBjb2RlX2Rvd25sb2FkOiB5ZXMKICAgIHNtb290aF9zY3JvbGw6IHllcwogICAgdGhlbWU6IGx1bWVuCiAgcGRmX2RvY3VtZW50OiAKICAgIHRvYzogeWVzCiAgICB0b2NfZGVwdGg6IDQKICAgIGZpZ19jYXB0aW9uOiB5ZXMKICAgIG51bWJlcl9zZWN0aW9uczogeWVzCiAgICBmaWdfd2lkdGg6IDMKICAgIGZpZ19oZWlnaHQ6IDMKICB3b3JkX2RvY3VtZW50OiAKICAgIHRvYzogeWVzCiAgICB0b2NfZGVwdGg6IDQKICAgIGZpZ19jYXB0aW9uOiB5ZXMKICAgIGtlZXBfbWQ6IHllcwplZGl0b3Jfb3B0aW9uczogCiAgY2h1bmtfb3V0cHV0X3R5cGU6IGlubGluZQotLS0KCmBgYHtjc3MsIGVjaG8gPSBGQUxTRX0KI1RPQzo6YmVmb3JlIHsKICBjb250ZW50OiAiVGFibGUgb2YgQ29udGVudHMiOwogIGZvbnQtd2VpZ2h0OiBib2xkOwogIGZvbnQtc2l6ZTogMS4yZW07CiAgZGlzcGxheTogYmxvY2s7CiAgY29sb3I6IG5hdnk7CiAgbWFyZ2luLWJvdHRvbTogMTBweDsKfQoKCmRpdiNUT0MgbGkgeyAgICAgLyogdGFibGUgb2YgY29udGVudCAgKi8KICAgIGxpc3Qtc3R5bGU6dXBwZXItcm9tYW47CiAgICBiYWNrZ3JvdW5kLWltYWdlOm5vbmU7CiAgICBiYWNrZ3JvdW5kLXJlcGVhdDpub25lOwogICAgYmFja2dyb3VuZC1wb3NpdGlvbjowOwp9CgpoMS50aXRsZSB7ICAgIC8qIGxldmVsIDEgaGVhZGVyIG9mIHRpdGxlICAqLwogIGZvbnQtc2l6ZTogMjJweDsKICBmb250LXdlaWdodDogYm9sZDsKICBjb2xvcjogRGFya1JlZDsKICB0ZXh0LWFsaWduOiBjZW50ZXI7CiAgZm9udC1mYW1pbHk6ICJHaWxsIFNhbnMiLCBzYW5zLXNlcmlmOwp9CgpoNC5hdXRob3IgeyAvKiBIZWFkZXIgNCAtIGFuZCB0aGUgYXV0aG9yIGFuZCBkYXRhIGhlYWRlcnMgdXNlIHRoaXMgdG9vICAqLwogIGZvbnQtc2l6ZTogMTVweDsKICBmb250LXdlaWdodDogYm9sZDsKICBmb250LWZhbWlseTogc3lzdGVtLXVpOwogIGNvbG9yOiBuYXZ5OwogIHRleHQtYWxpZ246IGNlbnRlcjsKfQoKaDQuZGF0ZSB7IC8qIEhlYWRlciA0IC0gYW5kIHRoZSBhdXRob3IgYW5kIGRhdGEgaGVhZGVycyB1c2UgdGhpcyB0b28gICovCiAgZm9udC1zaXplOiAxOHB4OwogIGZvbnQtd2VpZ2h0OiBib2xkOwogIGZvbnQtZmFtaWx5OiAiR2lsbCBTYW5zIiwgc2Fucy1zZXJpZjsKICBjb2xvcjogRGFya0JsdWU7CiAgdGV4dC1hbGlnbjogY2VudGVyOwp9CgpoMSB7IC8qIEhlYWRlciAxIC0gYW5kIHRoZSBhdXRob3IgYW5kIGRhdGEgaGVhZGVycyB1c2UgdGhpcyB0b28gICovCiAgICBmb250LXNpemU6IDIwcHg7CiAgICBmb250LXdlaWdodDogYm9sZDsKICAgIGZvbnQtZmFtaWx5OiAiVGltZXMgTmV3IFJvbWFuIiwgVGltZXMsIHNlcmlmOwogICAgY29sb3I6IGRhcmtyZWQ7CiAgICB0ZXh0LWFsaWduOiBjZW50ZXI7Cn0KCmgyIHsgLyogSGVhZGVyIDIgLSBhbmQgdGhlIGF1dGhvciBhbmQgZGF0YSBoZWFkZXJzIHVzZSB0aGlzIHRvbyAgKi8KICAgIGZvbnQtc2l6ZTogMThweDsKICAgIGZvbnQtd2VpZ2h0OiBib2xkOwogICAgZm9udC1mYW1pbHk6ICJUaW1lcyBOZXcgUm9tYW4iLCBUaW1lcywgc2VyaWY7CiAgICBjb2xvcjogbmF2eTsKICAgIHRleHQtYWxpZ246IGxlZnQ7Cn0KCmgzIHsgLyogSGVhZGVyIDMgLSBhbmQgdGhlIGF1dGhvciBhbmQgZGF0YSBoZWFkZXJzIHVzZSB0aGlzIHRvbyAgKi8KICAgIGZvbnQtc2l6ZTogMTZweDsKICAgIGZvbnQtd2VpZ2h0OiBib2xkOwogICAgZm9udC1mYW1pbHk6ICJUaW1lcyBOZXcgUm9tYW4iLCBUaW1lcywgc2VyaWY7CiAgICBjb2xvcjogbmF2eTsKICAgIHRleHQtYWxpZ246IGxlZnQ7Cn0KCmg0IHsgLyogSGVhZGVyIDQgLSBhbmQgdGhlIGF1dGhvciBhbmQgZGF0YSBoZWFkZXJzIHVzZSB0aGlzIHRvbyAgKi8KICAgIGZvbnQtc2l6ZTogMTRweDsKICBmb250LXdlaWdodDogYm9sZDsKICAgIGZvbnQtZmFtaWx5OiAiVGltZXMgTmV3IFJvbWFuIiwgVGltZXMsIHNlcmlmOwogICAgY29sb3I6IGRhcmtyZWQ7CiAgICB0ZXh0LWFsaWduOiBsZWZ0Owp9CgovKiBBZGQgZG90cyBhZnRlciBudW1iZXJlZCBoZWFkZXJzICovCi5oZWFkZXItc2VjdGlvbi1udW1iZXI6OmFmdGVyIHsKICBjb250ZW50OiAiLiI7Cgpib2R5IHsgYmFja2dyb3VuZC1jb2xvcjp3aGl0ZTsgfQoKLmhpZ2hsaWdodG1lIHsgYmFja2dyb3VuZC1jb2xvcjp5ZWxsb3c7IH0KCnAgeyBiYWNrZ3JvdW5kLWNvbG9yOndoaXRlOyB9Cgp9CmBgYAoKYGBge3Igc2V0dXAsIGluY2x1ZGU9RkFMU0V9CiMgY29kZSBjaHVuayBzcGVjaWZpZXMgd2hldGhlciB0aGUgUiBjb2RlLCB3YXJuaW5ncywgYW5kIG91dHB1dCAKIyB3aWxsIGJlIGluY2x1ZGVkIGluIHRoZSBvdXRwdXQgZmlsZXMuCmlmICghcmVxdWlyZSgia25pdHIiKSkgewogICBpbnN0YWxsLnBhY2thZ2VzKCJrbml0ciIpCiAgIGxpYnJhcnkoa25pdHIpCn0KaWYgKCFyZXF1aXJlKCJwYW5kZXIiKSkgewogICBpbnN0YWxsLnBhY2thZ2VzKCJwYW5kZXIiKQogICBsaWJyYXJ5KHBhbmRlcikKfQppZiAoIXJlcXVpcmUoImdncGxvdDIiKSkgewogIGluc3RhbGwucGFja2FnZXMoImdncGxvdDIiKQogIGxpYnJhcnkoZ2dwbG90MikKfQppZiAoIXJlcXVpcmUoInRpZHl2ZXJzZSIpKSB7CiAgaW5zdGFsbC5wYWNrYWdlcygidGlkeXZlcnNlIikKICBsaWJyYXJ5KHRpZHl2ZXJzZSkKfQoKaWYgKCFyZXF1aXJlKCJwbG90bHkiKSkgewogIGluc3RhbGwucGFja2FnZXMoInBsb3RseSIpCiAgbGlicmFyeShwbG90bHkpCn0KIyMjIwprbml0cjo6b3B0c19jaHVuayRzZXQoZWNobyA9IFRSVUUsICAgICAgICMgaW5jbHVkZSBjb2RlIGNodW5rIGluIHRoZSBvdXRwdXQgZmlsZQogICAgICAgICAgICAgICAgICAgICAgd2FybmluZyA9IEZBTFNFLCAgICMgc29tZXRpbWVzLCB5b3UgY29kZSBtYXkgcHJvZHVjZSB3YXJuaW5nIG1lc3NhZ2VzLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICMgeW91IGNhbiBjaG9vc2UgdG8gaW5jbHVkZSB0aGUgd2FybmluZyBtZXNzYWdlcyBpbgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICMgdGhlIG91dHB1dCBmaWxlLiAKICAgICAgICAgICAgICAgICAgICAgIHJlc3VsdHMgPSBUUlVFLCAgICAjIHlvdSBjYW4gYWxzbyBkZWNpZGUgd2hldGhlciB0byBpbmNsdWRlIHRoZSBvdXRwdXQKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAjIGluIHRoZSBvdXRwdXQgZmlsZS4KICAgICAgICAgICAgICAgICAgICAgIG1lc3NhZ2UgPSBGQUxTRSwKICAgICAgICAgICAgICAgICAgICAgIGNvbW1lbnQgPSBOQQogICAgICAgICAgICAgICAgICAgICAgKSAgCmBgYAogCiBcCiAKIyMgKipBc3NpZ25tZW50IE9iamVjdGl2ZXMqKiAKCiogTWFzdGVyIHRoZSBmdW5kYW1lbnRhbCBjb25jZXB0cyBvZiBwb2ludCBlc3RpbWF0aW9uIGFuZCBwZXJmb3JtYW5jZSBtZXRyaWNzCgoqIFVuZGVyc3RhbmQgdGhlIHRoZW9yZXRpY2FsIGZvdW5kYXRpb24gb2YgdGhlIG1ldGhvZCBvZiBtb21lbnRzIGVzdGltYXRvciAoTU1FKQoKKiBJbXBsZW1lbnQgTU1FIGluIFIsIGluY29ycG9yYXRpbmcgbnVtZXJpY2FsIGFwcHJveGltYXRpb24gbWV0aG9kcwoKXAoKKipMb2ctbG9naXN0aWMgRGlzdHJpYnV0aW9uKioKClRoZSBsb2ctbG9naXN0aWMgZGlzdHJpYnV0aW9uIChhbHNvIGtub3duIGFzIHRoZSBGaXNrIGRpc3RyaWJ1dGlvbikgaXMgYSBjb250aW51b3VzIHByb2JhYmlsaXR5IGRpc3RyaWJ1dGlvbiB0aGF0IGlzIHBhcnRpY3VsYXJseSB1c2VmdWwgaW4gY29udGV4dHMgd2hlcmUgZGF0YSBleGhpYml0IG5vbi1uZWdhdGl2ZSwgc2tld2VkIGJlaGF2aW9yIGFuZCB3aGVyZSB0aGUgaGF6YXJkIHJhdGUgaXMgdW5pbW9kYWwgKGluY3JlYXNlcyB0byBhIHBlYWsgYW5kIHRoZW4gZGVjcmVhc2VzKS4gSXQgaGFzIGJlZW4gd2lkZWx5IHVzZWQgaW4gdGhlIGFyZWFzIHN1Y2ggYXMgc3Vydml2YWwgYW5hbHlzaXMgYW5kIHJlbGlhYmlsaXR5IGVuZ2luZWVyaW5nLCBlbnZpcm9ubWVudGFsIHNjaWVuY2UsIGVjb25vbWljcywgcGhhcm1hY29sb2d5LCBmaW5hbmNlIGFuZCByaXNrIG1hbmFnZW1lbnQsIGV0Yy4gCgpGb3IgZ2l2ZW4gc2hhcGUgcGFyYW1ldGVyICRcYmV0YSQgYW5kIHNjYWxlIHBhcmFtZXRlciAkXGFscGhhJCwgdGhlIGN1bXVsYXRpdmUgZGlzdHJpYnV0aW9uIGZ1bmN0aW9uCgokJApGKHgpID0gXGZyYWN7MX17MSsoeC9cYWxwaGEpXnstXGJldGF9fQokJAoKQXMgYW4gZXhlcmNpc2UsIHlvdSBjYW4gZGVyaXZlIHRoZSBkZW5zaXR5IGluIHRoZSBmb2xsb3dpbmcgZm9ybQoKJCQKZih4KSA9IFxmcmFjeyhcYmV0YS9cYWxwaGEpKHgvXGFscGhhKV57XGJldGEtMX19e1sxKyh4L1xhbHBoYSleXGJldGFdXjJ9LCBcIFwgXHRleHR7IGZvciB9IFwgXCB4ID4gMC4KJCQKCkFmdGVyIHNvbWUgYWxnZWJyYSwgd2UgY2FuIGZpbmQgdGhlICRrJHRoIG1vbWVudAoKJCQKXG11X2sgPSBFW1hea10gPSBcYWxwaGFeayBCXGxlZnQoMStcZnJhY3trfXtcYmV0YX0sIDEgLSBcZnJhY3trfXtcYmV0YX0gXHJpZ2h0KS4KJCQKClRoaXMgYXNzaWdubWVudCB3aWxsIGZvY3VzIG9uIGZpbmRpbmcgTU1FIG9mIHBhcmFtZXRlcnMgJFxhbHBoYSQgYW5kICRcYmV0YSQgYmFzZWQgb24gYSByZWFsLXdvcmxkIGFwcGxpY2F0aW9uIGRhdGEgc2V0LgoKClwKCiMjICoqUXVlc3Rpb24gMTogRGVyaXZlIHRoZSBsb2ctbG9naXN0aWMgZGVuc2l0eSBmdW5jdGlvbiAqKgoKR2l2ZW4gdGhlIENERiBvZiB0aGUgdHdvLXBhcmFtZXRlciBsb2ctbG9naXN0aWMgZGlzdHJpYnV0aW9uCgokJApGKHgpID0gXGZyYWN7MX17MSsoeC9cYWxwaGEpXnstXGJldGF9fS4KJCQKJCQKXGJlZ2lue2FsaWduZWR9CkYoeCkKJj0gXGZyYWN7MX17MSArICh4L1xhbHBoYSleey1cYmV0YX19IFxcWzZwdF0KJj0gXGxlZnQoMSArICh4L1xhbHBoYSleey1cYmV0YX1ccmlnaHQpXnstMX0gXFxbMTBwdF0KCmYoeCkKJj0gXGZyYWN7ZH17ZHh9Rih4KSBcXFs2cHRdCiY9IC1cbGVmdCgxICsgKHgvXGFscGhhKV57LVxiZXRhfVxyaWdodCleey0yfQogICAgXGNkb3QKICAgIFxmcmFje2R9e2R4fVxsZWZ0KDEgKyAoeC9cYWxwaGEpXnstXGJldGF9XHJpZ2h0KSBcXFsxMHB0XQoKJj0gLVxsZWZ0KDEgKyAoeC9cYWxwaGEpXnstXGJldGF9XHJpZ2h0KV57LTJ9CiAgICBcY2RvdAogICAgXGxlZnQoMCArICgtXGJldGEpKHgvXGFscGhhKV57LVxiZXRhLTF9XGNkb3QgXGZyYWN7MX17XGFscGhhfVxyaWdodCkgXFxbMTBwdF0KCiY9IFxmcmFje1xiZXRhICh4L1xhbHBoYSleey0oXGJldGErMSl9fQogICAgICAgIHtcYWxwaGEgXGxlZnQoMSArICh4L1xhbHBoYSleey1cYmV0YX1ccmlnaHQpXjJ9IFxcWzEwcHRdCgomPSBcZnJhY3tcYmV0YSAoeC9cYWxwaGEpXntcYmV0YS0xfX0KICAgICAgICB7XGFscGhhIFxsZWZ0KDEgKyAoeC9cYWxwaGEpXntcYmV0YX1ccmlnaHQpXjJ9IFxcWzEycHRdCiY9IFxmcmFjeyhcYmV0YS9cYWxwaGEpKHgvXGFscGhhKV57XGJldGEtMX19e1sxKyh4L1xhbHBoYSleXGJldGFdXjJ9LCBcIFwgXHRleHR7IGZvciB9IFwgXCB4ID4gMC4KXGVuZHthbGlnbmVkfQokJApcCgojIyAqKlF1ZXN0aW9uIDI6IERpc3RyaWJ1dGlvbiBvZiBSZWNvdmVyeSBUaW1lIGZyb20gQSBTdXJnZXJ5KioKClRpbWUgdG8gcmVjb3ZlcnkgKGluIGRheXMpIGFmdGVyIGEgc3BlY2lmaWMga25lZSBzdXJnZXJ5IHByb2NlZHVyZS4gVGhpcyBmb2xsb3dzIGEgdHlwaWNhbCAqKmxvZy1sb2dpc3RpYyBwYXR0ZXJuKiogaW4gbWVkaWNhbCBzdXJ2aXZhbC9yZWNvdmVyeSBhbmFseXNpczoKCmBgYAo4LjIzLCAxMi43NCwgMTQuODMsIDE2LjYxLCAxOC4xNiwgMTkuNTUsIDIwLjgwLCAyMS45NCwgMjMuMDAsIDIzLjk4LCAyNC44OSwgMjUuNzUsIDI2LjU2LCAKMjcuMzQsIDI4LjA4LCAyOC43OSwgMjkuNDgsIDMwLjE1LCAzMC44MSwgMzEuNDUsIDMyLjA4LCAzMi43MCwgMzMuMzEsIDMzLjkyLCAzNC41MywgMzUuMTMsIAozNS43MywgMzYuMzMsIDM2LjkzLCAzNy41MywgMzguMTQsIDM4Ljc1LCAzOS4zNywgNDAuMDAsIDQwLjY0LCA0MS4yOSwgNDEuOTUsIDQyLjYzLCA0My4zMywgCjQ0LjA1LCA0NC43OSwgNDUuNTYsIDQ2LjM2LCA0Ny4yMCwgNDguMDgsIDQ5LjAyLCA1MC4wMywgNTEuMTIsIDUyLjMyLCA1My42NQpgYGAKQmFzZWQgb24gdGhlIGFib3ZlIGRhdGEgdG8gcGVyZm9ybSB0aGUgZm9sbG93aW5nIGFuYWx5c2lzLgoKYSkgVXNpbmcgbWV0aG9kIG9mIG1vbWVudCBlc3RpbWF0aW9uIHRvIGVzdGltYXRlICRcYWxwaGEkIGFuZCAkXGJldGEkLCBkZW5vdGVkIGJ5ICRcaGF0e1xhbHBoYX0kIGFuZCAkXGhhdHtcYmV0YX0kLCByZXNwZWN0aXZlbHkuCgpNb21lbnQgZ2VuZXJhdGluZyBmdW5jdGlvbiBmb3IgYSBsb2cgbG9naXN0aWMgZGlzdHJpYnV0aW9uLgokJApcbXVfayA9IEVbWF5rXSA9IFxhbHBoYV5rIEJcbGVmdCgxK1xmcmFje2t9e1xiZXRhfSxcLCAxLVxmcmFje2t9e1xiZXRhfVxyaWdodCkKJCQKVGhlIGZpcnN0IHR3byBtb21lbnRzIGFyZQokJApcbXVfMSA9IEVbWF0gPSBcYWxwaGEgQlxsZWZ0KDErXGZyYWN7MX17XGJldGF9LFwsIDEtXGZyYWN7MX17XGJldGF9XHJpZ2h0KQokJAokJApcbXVfMiA9IEVbWF4yXSA9IFxhbHBoYV4yIEJcbGVmdCgxK1xmcmFjezJ9e1xiZXRhfSxcLCAxLVxmcmFjezJ9e1xiZXRhfVxyaWdodCkKJCQKQmV0YSBGdW5jdGlvbiBJZGVudGl0eQokJApCKGEsYik9XGZyYWN7XEdhbW1hKGEpXEdhbW1hKGIpfXtcR2FtbWEoYStiKX0KJCQKU2ltcGxpZnlpbmcgYm90aCBtb21lbnRzCiQkClxtdV8xCj0KXGFscGhhClxHYW1tYVxsZWZ0KDErXGZyYWN7MX17XGJldGF9XHJpZ2h0KQpcR2FtbWFcbGVmdCgxLVxmcmFjezF9e1xiZXRhfVxyaWdodCkKJCQKJCQKXG11XzIKPQpcYWxwaGFeMgpcR2FtbWFcbGVmdCgxK1xmcmFjezJ9e1xiZXRhfVxyaWdodCkKXEdhbW1hXGxlZnQoMS1cZnJhY3syfXtcYmV0YX1ccmlnaHQpCiQkClRoZSBmaXJzdCB0d28gc2FtcGxlIG1vbWVudHMgYXJlCiQkCm1fMSA9IFxmcmFjezF9e259XHN1bV97aT0xfV5uIFhfaSA9IDM0LjE5MjIKJCQKJCQKbV8yID0gXGZyYWN7MX17bn1cc3VtX3tpPTF9Xm4gWF9pXjIgPSAxMjg4Ljg0NQokJApTdWJzdGl0dXRlIHRoZW9yZXRpY2FsIHNhbXBsZSBtb21lbnRzIGZvciBtb21lbnRzIHNpbmNlIHRoZXkncmUgZXF1YWwKJCQKMzQuMTkyMgo9ClxhbHBoYQpcR2FtbWFcbGVmdCgxK1xmcmFjezF9e1xiZXRhfVxyaWdodCkKXEdhbW1hXGxlZnQoMS1cZnJhY3sxfXtcYmV0YX1ccmlnaHQpCiQkCgokJAoxMjg4Ljg0NQo9ClxhbHBoYV4yClxHYW1tYVxsZWZ0KDErXGZyYWN7Mn17XGJldGF9XHJpZ2h0KQpcR2FtbWFcbGVmdCgxLVxmcmFjezJ9e1xiZXRhfVxyaWdodCkKJCQKU29sdmUgZmlyc3QgZXF1YXRpb24gZm9yICRcYWxwaGEkCiQkClxhbHBoYQo9ClxmcmFjezM0LjE5MjJ9CntcR2FtbWFcbGVmdCgxK1xmcmFjezF9e1xiZXRhfVxyaWdodCkKXEdhbW1hXGxlZnQoMS1cZnJhY3sxfXtcYmV0YX1ccmlnaHQpfQokJAoKUGx1ZyAkXGFscGhhJCBpbnRvIHNlY29uZCBlcXVhdGlvbgokJAoxMjg4Ljg0NQo9ClxsZWZ0KApcZnJhY3szNC4xOTIyfQp7XEdhbW1hXGxlZnQoMStcZnJhY3sxfXtcYmV0YX1ccmlnaHQpClxHYW1tYVxsZWZ0KDEtXGZyYWN7MX17XGJldGF9XHJpZ2h0KX0KXHJpZ2h0KV4yClxHYW1tYVxsZWZ0KDErXGZyYWN7Mn17XGJldGF9XHJpZ2h0KQpcR2FtbWFcbGVmdCgxLVxmcmFjezJ9e1xiZXRhfVxyaWdodCkKJCQKCkRpdmlkZSBib3RoIHNpZGVzIGJ5ICR7bV8xXjJ9JCAocmVwbGFjaW5nIHRoZW9yZXRpY2FsIHNhbXBsZSBtb21lbnRzIHdpdGggJHttXzF9JCBhbmQgJHttXzJ9JCkKJCQKXGZyYWN7bV8yfXttXzFeMn0KPQpcZnJhY3sKXEdhbW1hXGxlZnQoMStcZnJhY3syfXtcYmV0YX1ccmlnaHQpClxHYW1tYVxsZWZ0KDEtXGZyYWN7Mn17XGJldGF9XHJpZ2h0KQp9ewpcbGVmdFsKXEdhbW1hXGxlZnQoMStcZnJhY3sxfXtcYmV0YX1ccmlnaHQpClxHYW1tYVxsZWZ0KDEtXGZyYWN7MX17XGJldGF9XHJpZ2h0KQpccmlnaHRdXjIKfQokJAoKRGVmaW5lCiQkCmcoXGJldGEpCj0KXGZyYWN7ClxHYW1tYVxsZWZ0KDErXGZyYWN7Mn17XGJldGF9XHJpZ2h0KQpcR2FtbWFcbGVmdCgxLVxmcmFjezJ9e1xiZXRhfVxyaWdodCkKfXsKXGxlZnRbClxHYW1tYVxsZWZ0KDErXGZyYWN7MX17XGJldGF9XHJpZ2h0KQpcR2FtbWFcbGVmdCgxLVxmcmFjezF9e1xiZXRhfVxyaWdodCkKXHJpZ2h0XV4yCn0KLQpcZnJhY3ttXzJ9e21fMV4yfQokJApTb2x2ZSAkZyhcYmV0YSk9MCQgdG8gb2J0YWluICRcaGF0e1xiZXRhfSQsIGFuZCB0aGVuIHN1YnN0aXR1dGUgaW50bwokJApcaGF0e1xhbHBoYX0KPQpcZnJhY3ttXzF9CntcR2FtbWFcbGVmdCgxK1xmcmFjezF9e1xoYXR7XGJldGF9fVxyaWdodCkKXEdhbW1hXGxlZnQoMS1cZnJhY3sxfXtcaGF0e1xiZXRhfX1ccmlnaHQpfQokJAoKQ29kZSB0byBzb2x2ZSAkZyhcYmV0YSk9MCQgYW5kIG9idGFpbiBib3RoICRcaGF0e1xiZXRhfSQgYW5kICRcaGF0e1xhbHBoYX0kLgpgYGB7cn0KeCA8LSBjKDguMjMsIDEyLjc0LCAxNC44MywgMTYuNjEsIDE4LjE2LCAxOS41NSwgMjAuODAsIDIxLjk0LCAyMy4wMCwgMjMuOTgsIDI0Ljg5LCAyNS43NSwgMjYuNTYsIAoyNy4zNCwgMjguMDgsIDI4Ljc5LCAyOS40OCwgMzAuMTUsIDMwLjgxLCAzMS40NSwgMzIuMDgsIDMyLjcwLCAzMy4zMSwgMzMuOTIsIDM0LjUzLCAzNS4xMywgCjM1LjczLCAzNi4zMywgMzYuOTMsIDM3LjUzLCAzOC4xNCwgMzguNzUsIDM5LjM3LCA0MC4wMCwgNDAuNjQsIDQxLjI5LCA0MS45NSwgNDIuNjMsIDQzLjMzLCAKNDQuMDUsIDQ0Ljc5LCA0NS41NiwgNDYuMzYsIDQ3LjIwLCA0OC4wOCwgNDkuMDIsIDUwLjAzLCA1MS4xMiwgNTIuMzIsIDUzLjY1KQoKIzFzdCBhbmQgMm5kIG1vbWVudHMKbTEgPSBtZWFuKHgpCm0yID0gbWVhbih4XjIpCgojcm9vdCBmaW5kaW5nIGVxdWF0aW9uCmdiIDwtIGZ1bmN0aW9uKEIpewogICggKGdhbW1hKDEgKyAyL0IpICogZ2FtbWEoMSAtIDIvQikpIC8gKChnYW1tYSgxICsgMS9CKSAqIGdhbW1hKDEgLSAxL0IpKV4yKSApIC0gKG0yIC8gbTFeMikKfQoKIyBzaW5jZSBiZXRhID4gMgp4X3ZhbHMgPC0gc2VxKDIuMDEsIDEwLCBsZW5ndGgub3V0ID0gMjAwKQpkZnJhbWUgPC0gZGF0YS5mcmFtZSh4ID0geF92YWxzLCB5ID0gZ2IoeF92YWxzKSkKCmdrLnBsb3QgPC0gZ2dwbG90KGRmcmFtZSwgYWVzKHggPSB4LCB5ID0geSkpICsKICBnZW9tX2xpbmUoY29sb3IgPSAic3RlZWxibHVlIiwgc2l6ZSA9IDEpICsKICBnZW9tX2hsaW5lKHlpbnRlcmNlcHQgPSAwLCBsaW5ldHlwZSA9ICJkYXNoZWQiLCBhbHBoYSA9IDAuNSkgKyAgIyB4LWF4aXMKICBnZW9tX3ZsaW5lKHhpbnRlcmNlcHQgPSAwLCBsaW5ldHlwZSA9ICJkYXNoZWQiLCBhbHBoYSA9IDAuNSkgKyAgIyB5LWF4aXMKICBsYWJzKHRpdGxlID0gIlRoZSBjdXJ2ZSBvZiBmdW5jdGlvbiBnKEJldGEpIiwKICAgICAgIHggPSAiQmV0YSIsIAogICAgICAgeSA9ICJnKEJldGEpIikgKwogICAgdGhlbWUocGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChoanVzdCA9IDAuNSksCiAgICAgICAgcGxvdC5tYXJnaW4gPSBtYXJnaW4odCA9IDM1LCByID0gMjAsIGIgPSAzMCwgbCA9IDMwLCB1bml0ID0gInB0IikpCmdncGxvdGx5KGdrLnBsb3QpCmBgYAoKYGBge3J9CiNmaW5kIHRoZSByb290IChiZXRhIHZhbHVlIHdoZW4gZXF1YXRpb24gPSAwKQpCIDwtIHVuaXJvb3QoZ2IsIGludGVydmFsID0gYyg1LCA3KSkkcm9vdAoKI3VzZSB0aGUgYmV0YSB0byBmaW5kIHRoZSBhbHBoYQphbHBoYSA8LSBtMS8gKGdhbW1hKDErMS9CKSAqIGdhbW1hKDEtMS9CKSkKcGFuZGVyKGNiaW5kKGJldGEgPSBCLCBhbHBoYSA9IGFscGhhKSkKYGBgCgpiKSBTaW5jZSB0aGUgbW9tZW50IGVzdGltYXRlcyAkXGhhdHtcYWxwaGF9JCBhbmQgJFxoYXR7XGJldGF9JCBhcmUgcmFuZG9tLCBjb25zdHJ1Y3QgYm9vdHN0cmFwIHNhbXBsaW5nIGRpc3RyaWJ1dGlvbnMgZm9yIGVhY2guIFRvIHZpc3VhbGl6ZSB0aGVzZSBkaXN0cmlidXRpb25zLCBwbG90IHNlcGFyYXRlIGJvb3RzdHJhcCBoaXN0b2dyYW1zIGZvciAkXGhhdHtcYWxwaGF9JCBhbmQgJFxoYXR7XGJldGF9JC4gIGhlbiwgb3ZlcmxheSBhIHNtb290aCBkZW5zaXR5IGN1cnZlIG9uIGVhY2ggaGlzdG9ncmFtIHVzaW5nIEdhdXNzaWFuIGtlcm5lbCBkZW5zaXR5IGVzdGltYXRpb24uIEZpbmFsbHksIGRlc2NyaWJlIHRoZSBwYXR0ZXJucyBvZiB0aGVzZSBkZW5zaXR5IGN1cnZlcy4KCmBgYHtyfQpib290c3RyYXAgPC0gZnVuY3Rpb24gKGRhdGEsIHN0YXRpc3RpY3MsIEIpewogIG4gPC0gbGVuZ3RoKGRhdGEpCiAgc3RhdHMgPC0gbWF0cml4KG5yb3cgPSBCLCBuY29sID0gMikKICBmb3IgKGkgaW4gMTpCKXsKICAgIG5ld19kYXRhIDwtIHNhbXBsZShkYXRhLCBuLCByZXBsYWNlID0gVFJVRSkKICAgIHN0YXRzW2ksXSA8LSBzdGF0aXN0aWNzKG5ld19kYXRhKQogIH0KICByZXR1cm4oc3RhdHMpCn0KICAKc3RhdGlzdGljcyA8LSBmdW5jdGlvbiAoZGF0YSl7CiAgbTEgPSBtZWFuKGRhdGEpCiAgbTIgPSBtZWFuKGRhdGFeMikKICAKICBnYiA8LSBmdW5jdGlvbihCKXsKICAoIChnYW1tYSgxICsgMi9CKSAqIGdhbW1hKDEgLSAyL0IpKSAvICgoZ2FtbWEoMSArIDEvQikgKiBnYW1tYSgxIC0gMS9CKSleMikgKSAtIChtMiAvIG0xXjIpfQogIAogIHhfdmFscyA8LSBzZXEoMi4wMSwgMjAwLCBsZW5ndGgub3V0ID0gMTAwKQogIGRmcmFtZSA8LSBkYXRhLmZyYW1lKHggPSB4X3ZhbHMsIHkgPSBnYih4X3ZhbHMpKQogIAogIEIgPC0gdW5pcm9vdChnYiwgaW50ZXJ2YWwgPSBjKDMsIDE1KSkkcm9vdAogIGFscGhhIDwtIG0xLyAoZ2FtbWEoMSsxL0IpICogZ2FtbWEoMS0xL0IpKQogIAogIGMoQiwgYWxwaGEpCn0KCnN0YXRzIDwtIGJvb3RzdHJhcCh4LCBzdGF0aXN0aWNzLCA1MDAwKQoKYmV0YV9ib290ICA8LSBzdGF0c1ssIDFdCmFscGhhX2Jvb3QgPC0gc3RhdHNbLCAyXQoKaGlzdChiZXRhX2Jvb3QsCiAgICAgYnJlYWtzID0gMjAsCiAgICAgcHJvYmFiaWxpdHkgPSBUUlVFLAogICAgIHhsYWIgPSBleHByZXNzaW9uKGhhdChiZXRhKSksCiAgICAgbWFpbiA9ICJCb290c3RyYXAgU2FtcGxpbmcgRGlzdHJpYnV0aW9uXG4gb2YgQmV0YSBIYXQiLAogICAgIGNleC5tYWluID0gMC45LAogICAgIGNvbCA9ICJsaWdodGdyYXkiLAogICAgIGJvcmRlciA9ICJ3aGl0ZSIpCgpsaW5lcyhkZW5zaXR5KGJldGFfYm9vdCksIGNvbCA9ICJuYXZ5IiwgbHdkID0gMikKYGBgCgpUaGUgYm9vdHN0cmFwIHNhbXBsaW5nIGRpc3RyaWJ1dGlvbiBvZiAkXGhhdHtcYmV0YX0kIGlzIHJvdWdobHkgYmVsbC1zaGFwZWQgd2l0aCBhIGxpdHRsZSBiaXQgb2Ygc2tld25lc3MgdG8gdGhlIHJpZ2h0LiBJdOKAmXMgY2VudGVyZWQgY2xvc2UgdG8gdGhlIG9yaWdpbmFsIE1vTSBlc3RpbWF0ZSwgd2hpY2ggc3VnZ2VzdHMgdGhlIGVzdGltYXRvciBpcyBub3QgaGVhdmlseSBiaWFzZWQuIFRoZSBzcHJlYWQgb2YgdGhlIGhpc3RvZ3JhbSBzaG93cyBob3cgbXVjaCAkXGhhdHtcYmV0YX0kIHdvdWxkIHZhcnkgaWYgd2UgY29sbGVjdGVkIG5ldyBzYW1wbGVzIGZyb20gdGhlIHBvcHVsYXRpb24uIFNpbmNlIHRoZSBkaXN0cmlidXRpb24gaXNu4oCZdCBleHRyZW1lbHkgd2lkZSBvciBpcnJlZ3VsYXIsIGl0IHN1Z2dlc3RzIHRoYXQgdGhlIE1vTSBlc3RpbWF0b3IgZm9yICR7XGJldGF9JCBpcyBmYWlybHkgc3RhYmxlIGZvciB0aGlzIHNhbXBsZS4KYGBge3J9Cmhpc3QoYWxwaGFfYm9vdCwKICAgICBicmVha3MgPSAyMCwKICAgICBwcm9iYWJpbGl0eSA9IFRSVUUsCiAgICAgeGxhYiA9IGV4cHJlc3Npb24oaGF0KGFscGhhKSksCiAgICAgbWFpbiA9ICJCb290c3RyYXAgU2FtcGxpbmcgRGlzdHJpYnV0aW9uXG4gb2YgQWxwaGEgSGF0IiwKICAgICBjZXgubWFpbiA9IDAuOSwKICAgICBjb2wgPSAibGlnaHRncmF5IiwKICAgICBib3JkZXIgPSAid2hpdGUiKQoKbGluZXMoZGVuc2l0eShhbHBoYV9ib290KSwgY29sID0gImRhcmtyZWQiLCBsd2QgPSAyKQpgYGAKClRoZSBib290c3RyYXAgc2FtcGxpbmcgZGlzdHJpYnV0aW9uIG9mICRcaGF0e1xhbHBoYX0kIGlzIGFwcHJveGltYXRlbHkgc3ltbWV0cmljIGFuZCBiZWxsLXNoYXBlZC4gSXTigJlzIGNlbnRlcmVkIHByZXR0eSBjbG9zZSB0byB0aGUgb3JpZ2luYWwgTW9NIGVzdGltYXRlLCB3aGljaCBzdWdnZXN0cyB0aGUgZXN0aW1hdG9yIGlzIG5vdCBoZWF2aWx5IGJpYXNlZC4gVGhlIHNwcmVhZCBvZiB0aGUgaGlzdG9ncmFtIHNob3dzIGhvdyBtdWNoICRcaGF0e1xhbHBoYX0kIHdvdWxkIHZhcnkgaWYgd2UgY29sbGVjdGVkIG5ldyBzYW1wbGVzIGZyb20gdGhlIHBvcHVsYXRpb24uIFNpbmNlIHRoZSBkaXN0cmlidXRpb24gaXNu4oCZdCBleHRyZW1lbHkgd2lkZSBvciBpcnJlZ3VsYXIsIGl0IHN1Z2dlc3RzIHRoYXQgdGhlIE1vTSBlc3RpbWF0b3IgZm9yICR7XGFscGhhfSQgaXMgZmFpcmx5IHN0YWJsZSBmb3IgdGhpcyBzYW1wbGUuCg==